Opis analizy

Analizie poddany został zbiór danych dotyczących połowu śledzia atlantyckiego oraz warunków środowiska w jakim żyje. Celem analizy była próba odpowiedzi na pytanie dlaczego od pewnego czasu złowione ryby mają coraz mniejszy rozmiar. Eksploracja rozpoczyna się od załadowania odpowiednich bibliotek, następnie ustawione zostaje ziarno zapewniające powtarzalność uzyskanych wyników oraz wczytany zostaje zbiór danych. Sprawdzona zostaje liczba pustych wartości oraz kształt zbioru. Ze względu na relatywnie niską ilość poszczególnych wartości NA zostaje podjęta decyzja o ich usunięciu.

Później przyglądano się bliżej poszczególnym atrybutom, ich rozkładom wartości oraz postawowym statyskom (mean, min, max, sd, etc.).

Kolejnym krokiem było sprawdzenie wzajemnych zależności pomiędzy atrybutami (włączając w to zmienną przewidywaną). Większość atrybutów nie wykazuje znaczącej wzajemnej zależności. Jednak kilka zmiennych, w szczególności dotyczących dostępności planktonu, jest silnie skorelowane. Dla zmiennych wykazujących wzajemne powiązania dokonano dodatkowej wizualizacji. Żaden z atrybutów nie jest bezpośrednio skorelowany z wartością zmiennej przewidywanej.

Następnie sporządzono wykres zmiany długości śledzia w czasie, jasno wskazujący na spadek średniego rozmiaru złowionych ryb.

Następny etap opierał się na zbudowaniu regresora, mającego za cel predykcję długości śledzia. Dokonano podziału zbioru na część treningową, walidacyjną oraz testową. W pierwszej kolejności przetestowano prosty model regresji liniowej i zmierzono jego jakość predykcji przy pomocy metryk RMSE oraz R2. Drugim wybranym modelem był XGBoost - wydajna biblioteka udostępniająca rozwiązania z dziedziny gradient boosting, oparte na drzewach decyzyjnych. Jest często używany w problemach analizy danych. Model wymagał dodatkowej konwersji danych na kompatybilny format. Następnie przy pomocy zbioru walidacyjnego dobrano parametry, tak, aby zmniejszyć overfitting. Użycie XGBoost pozwoliło na wyraźną poprawę wyników predykcji.

Finalnie, dla modelu o najlepszej skuteczności(XGBoost), przy pomocy metody LIME dokonano analizy ważności atrybutów oraz ich wpływu na estymowaną wartość w poszczególnych przypadkach.

Do najważniejszych cech należą odwiednio:

  • sst - 0.57,
  • recr - 0.14,
  • lcop1 - 0.06.

XGBoost dokonuje automatycznej analizy zależności zmiennych, na wykresie feature importance można zauważyć, że tylko jedna z silnie skorelowanych zmiennych ma wyższe znaczenie dla modelu.

Z dużego wpływu temperatury przy powierzchni wody (sst) na podejmowaną przez model decyzję można by wnioskować, że drobny wzrost temperatury (mniej więcej od podobnego okresu, gdy długość śledzia zaczęła spadać) mógł mieć wpływ na zmniejszenie rozmiaru poławianych śledzi.

Wstępne przetwarzanie

Wykorzystane biblioteki

library(reshape2)
library(plyr)
library(dplyr)
library(caret)
library(corrplot)
library(ggplot2)
library(plotly)
library(glmnet)
library(gbm)
library(xgboost)
library(lime)
library(nnet)

Powtarzalność wyników

set.seed(10)

Wczytywanie danych

Wartości puste oznaczone są znakiem “?”, przy wczywaniu zbioru korzystamy z parametru na.strings do zastąpienia ich wartością NA

data <- read.csv("sledzie.csv", header=TRUE, na.strings = "?")

Brakujące dane według kolumn

nmissing <- function(x) sum(is.na(x))

knitr::kable(colwise(nmissing)(data))
X length cfin1 cfin2 chel1 chel2 lcop1 lcop2 fbar recr cumf totaln sst sal xmonth nao
0 0 1581 1536 1555 1556 1653 1591 0 0 0 0 1584 0 0 0

Wymiary zbioru: (52582, 16)

Całkowita liczba wierszy: 52582

Liczba wierszy zawierających pustą wartość: 10094

Liczba wierszy po usunięciu pustych wartości: 42488

Przykładowe wartości ze zbioru

X length cfin1 cfin2 chel1 chel2 lcop1 lcop2 fbar recr cumf totaln sst sal xmonth nao
24923 24922 24.0 0.00000 0.70118 11.50000 5.67765 23.00000 9.17171 0.227 421391 0.1480941 730351.2 13.63160 35.50835 2 -2.86
18378 18377 25.5 0.00000 0.72727 2.59316 33.12121 5.35088 37.23232 0.547 148045 0.2714717 179955.4 13.39520 35.49992 5 -2.14
22182 22181 25.5 0.00000 0.70118 11.50000 5.67765 23.00000 9.17171 0.227 421391 0.1480941 730351.2 13.63160 35.50835 4 -2.86
34827 34826 25.0 3.14462 5.80145 4.26148 18.43839 7.67148 27.47359 0.254 474983 0.2073709 534157.2 13.89253 35.50455 5 -0.75
46082 46081 22.5 0.05455 0.23134 1.48629 11.67093 1.83515 17.80753 0.591 473462 0.3758273 306067.6 14.35600 35.52329 9 0.72

Analiza wartości atrybutów

W zbiorze znajdują się następujące kolumny(atrybuty):

  • length: długość złowionego śledzia [cm],
  • cfin1: dostępność planktonu [zagęszczenie Calanus finmarchicus gat. 1],
  • cfin2: dostępność planktonu [zagęszczenie Calanus finmarchicus gat. 2],
  • chel1: dostępność planktonu [zagęszczenie Calanus helgolandicus gat. 1],
  • chel2: dostępność planktonu [zagęszczenie Calanus helgolandicus gat. 2],
  • lcop1: dostępność planktonu [zagęszczenie widłonogów gat. 1],
  • lcop2: dostępność planktonu [zagęszczenie widłonogów gat. 2],
  • fbar: natężenie połowów w regionie [ułamek pozostawionego narybku],
  • recr: roczny narybek [liczba śledzi],
  • cumf: łączne roczne natężenie połowów w regionie [ułamek pozostawionego narybku],
  • totaln: łączna liczba ryb złowionych w ramach połowu [liczba śledzi],
  • sst: temperatura przy powierzchni wody [°C],
  • sal: poziom zasolenia wody [Knudsen ppt],
  • xmonth: miesiąc połowu [numer miesiąca],
  • nao: oscylacja północnoatlantycka [mb].

Długość złowionego śledzia length

Na histogramie długości złowionego śledzia wyraźnie rysuje się rozkład normalny. Odchylenie standardowe jest relatywnie niskie, równe jest 1,65 cm. 

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. Std.Dev. 
##    19.00    24.00    25.50    25.30    26.50    32.50     1.65

Dostępność planktonu cfin1

Dostępność planktonu cfin1 w okresie pomiarów była zdecydowanie niska, rozkład wartości jest w ponad 80% zdominowany przez wartość minimalną. Zdarzają się nieliczne wysokie wartości. 

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. Std.Dev. 
##   0.0000   0.0000   0.1111   0.4457   0.3333  37.6667   0.9800

Dostępność planktonu cfin2

Podobnie jak w przypadku cfin1, rozkład cfin2 jest zdominowany przez niskie wartości, jednak w znacznie mniejszym stopniu. W tym przypadku również występują wyższe wartości cechy, jednak tutaj układają się w widoczne grupy. 

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. Std.Dev. 
##   0.0000   0.2778   0.7012   2.0269   1.7936  19.3958   3.7200

Dostępność planktonu chel1

Analogicznie jak pozostałych przypadkach widoczna jest niska dostępność planktonu chel1. 

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. Std.Dev. 
##    0.000    2.469    5.750   10.016   11.500   75.000   14.320

Dostępność planktonu chel2

Dostępność chel2 w całym zbiorze jest dość zróżnicowana, w rozkładzie występuje dużo różnych wartości o podobnej częstości występowania. 

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. Std.Dev. 
##    5.238   13.427   21.435   21.197   27.193   57.706    9.990

Dostępność planktonu lcop1

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. Std.Dev. 
##   0.3074   2.5479   7.0000  12.8386  21.2315 115.5833  15.1000

Dostępność planktonu lcop2

Podobnie jak chel2, występowanie planktonu lcop2 jest dość zróżnicowane. Wartości obu planktonów wydają się być skorelowane.

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. Std.Dev. 
##    7.849   17.808   24.859   28.396   37.232   68.736   13.870

Natężenie połowów w regionie fbar

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. Std.Dev. 
##   0.0680   0.2270   0.3320   0.3306   0.4650   0.8490   0.1600

Roczny narybek recr

Zmienna recr przedstawia roczny narybek reprezentowany w postaci liczby złowionych ryb, wartości liczbowe osiągają dość wysokie wartości, niejednokrotnie przekraczające milion. 

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. Std.Dev. 
##   140515   360061   421391   519877   724151  1565890   270639

Łączne roczne natężenie połowów w regionie cumf

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. Std.Dev. 
##  0.06833  0.14809  0.23191  0.22987  0.29803  0.39801  0.09000

Łączna liczba ryb złowionych w ramach połowu totaln

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. Std.Dev. 
##   144137   306068   539558   515082   730351  1015595   221389

Temperatura przy powierzchni wody sst

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. Std.Dev. 
##    12.77    13.60    13.86    13.87    14.16    14.73     0.44

Poziom zasolenia wody sal

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. Std.Dev. 
##    35.40    35.51    35.51    35.51    35.52    35.61     0.04

Miesiąc połowu xmonth

Dane dotyczące miesiąca mają w przybliżeniu rozkład normalny, większość pomiarów pochodzi z okresu letnio-jesiennego. 

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. Std.Dev. 
##    1.000    5.000    8.000    7.252    9.000   12.000    2.760

Oscylacja północnoatlantycka nao

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. Std.Dev. 
## -4.89000 -1.90000  0.20000 -0.09642  1.63000  5.08000  2.25000

Zależności między zmiennymi

Z wykresu można wywnioskować, że wiekszość zmiennych jest skorelowana ze sobą jedynie w lekkim stopniu lub wcale. Jest jednak kilka wyjątków, przede wszystkim w przypadku dostępności niektórych gatunków planktonu. Występują dwie bardzo silne korelacje:

  • 0.96 między chel1, a lcop1,
  • 0.89 między chel2, a lcop2.

oraz kolejna para silnych korelacji dodatnich:

  • 0.82 między fbar, a cumf
  • 0.65 między cfin2, a lcop2

i jedna silna korelacja ujemna:

  • -0.71 między totaln, a cumf

Na wykresie widzimy bardzo silnie skorelowane zmienne lcop1 oraz chel1

W tym przypadku mamy do czynienia z negatywną korelacją zmiennych totaln i cumf. Siła korelacji jest mniejsza niż w poprzednich przypadkach, jednak wciąż można łatwo dostrzec malejący trend.

Zmiana długości śledzia w czasie

Zaprezentowany wykres, celem zwiększenia czytelności, jest przygotowany na bazie próbki 1000 pomiarów. Dzięki wygładzonej linii średniej można łatwo dostrzec jaki był kształt trendu w czasie. Interaktywny wykres pozwala dokładnie prześledzić zmianę wartości, na wykresie wyraźnie widać, że po około 16000. próbce średnia długość śledzia zaczyna spadać, dzieje się tak aż do końca badanego okresu.

Budowa regresora oraz pomiar jakości predykcji

Podział zbioru

train_p <- 0.7
valid_p <- 0.5

idx_train <- createDataPartition(y = data_clean$length,  p = train_p, list = FALSE)

train_set <- data_clean[,-1][ idx_train,]
test_set  <- data_clean[,-1][-idx_train,]

idx_valid <- createDataPartition(y = test_set$length,  p = valid_p, list = FALSE)

valid_set <- test_set[ idx_valid,]
test_set  <- test_set[-idx_valid,]

Rozmiary zbiorów: treningowego, walidacyjnego oraz testowego

dim(train_set)
## [1] 29744    15
dim(valid_set)
## [1] 6372   15
dim(test_set)
## [1] 6372   15

Regresja liniowa

simple_regr <- lm(length~. ,data = train_set)
simple_prediction <- predict(simple_regr, test_set[,-1])

Podsumowanie parametrów modelu regresji lm

summary(simple_regr)
## 
## Call:
## lm(formula = length ~ ., data = train_set)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.8848 -0.9335  0.0018  0.9126  6.9429 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  6.487e+01  8.458e+00   7.669 1.79e-14 ***
## cfin1        1.384e-01  9.372e-03  14.769  < 2e-16 ***
## cfin2        1.698e-02  5.224e-03   3.250  0.00115 ** 
## chel1       -1.023e-02  2.250e-03  -4.546 5.50e-06 ***
## chel2       -4.628e-03  3.316e-03  -1.396  0.16284    
## lcop1        2.181e-02  2.097e-03  10.401  < 2e-16 ***
## lcop2        1.129e-02  2.756e-03   4.094 4.25e-05 ***
## fbar         6.126e+00  1.132e-01  54.135  < 2e-16 ***
## recr        -3.905e-07  3.729e-08 -10.470  < 2e-16 ***
## cumf        -1.021e+01  2.112e-01 -48.323  < 2e-16 ***
## totaln      -5.909e-07  7.139e-08  -8.277  < 2e-16 ***
## sst         -1.274e+00  2.674e-02 -47.659  < 2e-16 ***
## sal         -6.086e-01  2.395e-01  -2.542  0.01104 *  
## xmonth       9.068e-03  2.863e-03   3.167  0.00154 ** 
## nao          3.519e-02  5.582e-03   6.305 2.92e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.358 on 29729 degrees of freedom
## Multiple R-squared:  0.3224, Adjusted R-squared:  0.3221 
## F-statistic:  1010 on 14 and 29729 DF,  p-value: < 2.2e-16

Wyniki

R2 0.3156754

RMSE 1.3586153

XGBoost

Przygotowanie zbioru do wersji kompatybilnej z XGBoost

xgb_train_set <- xgb.DMatrix(data = as.matrix(train_set[,-1]), label = train_set$length)

xgb_valid_set <- xgb.DMatrix(data = as.matrix(valid_set[,-1]), label = valid_set$length)

xgb_test_set <- xgb.DMatrix(data = as.matrix(test_set[,-1]), label = test_set$length)


watchlist <- list(train=xgb_train_set, valid=xgb_valid_set)


xg <- xgb.train(data = xgb_train_set, max.depth=32, nthread = 5, nrounds = 20, watchlist=watchlist)
## [1]  train-rmse:17.424788    valid-rmse:17.426558 
## [2]  train-rmse:12.230176    valid-rmse:12.232804 
## [3]  train-rmse:8.605362 valid-rmse:8.609152 
## [4]  train-rmse:6.084421 valid-rmse:6.088549 
## [5]  train-rmse:4.341674 valid-rmse:4.347968 
## [6]  train-rmse:3.150393 valid-rmse:3.159253 
## [7]  train-rmse:2.353688 valid-rmse:2.366677 
## [8]  train-rmse:1.837287 valid-rmse:1.854480 
## [9]  train-rmse:1.517758 valid-rmse:1.538785 
## [10] train-rmse:1.331451 valid-rmse:1.356127 
## [11] train-rmse:1.228172 valid-rmse:1.255594 
## [12] train-rmse:1.173326 valid-rmse:1.203040 
## [13] train-rmse:1.144919 valid-rmse:1.176587 
## [14] train-rmse:1.130379 valid-rmse:1.163402 
## [15] train-rmse:1.122979 valid-rmse:1.156788 
## [16] train-rmse:1.119246 valid-rmse:1.153730 
## [17] train-rmse:1.117372 valid-rmse:1.152286 
## [18] train-rmse:1.116414 valid-rmse:1.151726 
## [19] train-rmse:1.115926 valid-rmse:1.151528 
## [20] train-rmse:1.115675 valid-rmse:1.151489

Wyniki

R2 0.5224293

RMSE 1.1351347

Analiza ważności atrybutów modelu XGBoost

Przy pomocy pakietu Lime można zaobserwować wpływ poszczególnych atrybutów na wartość przewidywanej zmiennej dla kilku wybranych próbek. Analizując kilka wybranych próbek można łatwo zauważyć pewne prawidłowości:

  • wysoka temperatura sst > 14.2 bardzo negatywnie wpływa rozmiar śledzia,
  • wyższa dostępność plaktonu lcop1 > 21.23 ma lekki pozytywny wpływ,
  • okres letni wpływa pozytywnie na długość, a jesienno-zimowy lekko negatywnie,
  • niskie natężenie połowów w rejonie ma negatywny wpływ na długość.

W przypadku większej ilości przypadków można je zwizualizować na wykresie w postaci heatmap.

XGBoost oferuje wbudowane feature importance mierzone przy pomocy parametru gain. Pośród użytych atrybutów zdecydowanie dominuje sst, mając ponad połowę z całego zysku przy podejmowaniu decyzji. Następny jest atrybut recr mający ponaddwukrotnie większą ważność od swojego następnika. Pozostałe zmienne mają już mniejszy, ale bardziej równomierny wkład w decyzję.

importanceRaw <- xgb.importance(feature_names = colnames(test_set[,-1]), model = xg, data = xgb_test_set)
xgb.plot.importance(importance_matrix = importanceRaw)

Z dużego wpływu temperatury przy powierzchni wody (sst) na podejmowaną przez model decyzję można by wnioskować, że drobny wzrost temperatury (mniej więcej od podobnego okresu, gdy długość śledzia zaczęła spadać) mógł mieć wpływ na zmniejszenie rozmiaru poławianych śledzi.